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We give an introduction to some of the numerical aspects in quantum chaos. The classical dynam- 
ics of two-dimensional area-preserving maps on the torus is illustrated using the standard map 
and a perturbed cat map. The quantization of area-preserving maps given by their generating 
function is discussed and for the computation of the eigenvalues a computer program in Python is 
presented. We illustrate the eigenvalue distribution for two types of perturbed cat maps, one lead- 
ing to COE and the other to CUE statistics. For the eigenfunctions of quantum maps we study 
the distribution of the eigenvectors and compare them with the corresponding random matrix 
distributions. The Husimi representation allows for a direct comparison of the localization of the 
eigenstates in phase space with the corresponding classical structures. Examples for a perturbed 
q ■ cat map and the standard map with different parameters are shown. 

Billiard systems and the corresponding quantum billiards are another important class of sys- 
tems (which are also relevant to applications, for example in mesoscopic physics). We provide a 
detailed exposition of the boundary integral method, which is one important method to determine 
the eigenvalues and eigenfunctions of the Helmholtz equation. We discuss several methods to de- 
^ . termine the eigenvalues from the Fredholm equation and illustrate them for the stadium billiard. 
The occurrence of spurious solutions is discussed in detail and illustrated for the circular billiard, 
the stadium billiard, and the annular sector billiard. 

We emphasize the role of the normal derivative function to compute the normalization of 
eigenfunctions, momentum representations or autocorrelation functions in a very efficient and 
direct way. Some examples for these quantities are given and discussed. 
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1 Introduction 

In this text, which is an expanded version of lectures held at a summer school in Bologna in 
2001, we give an introduction to some of the numerical aspects in quantum chaos; some of the 
sections on the boundary integral method contain more advanced material. In quantum chaos one 
studies quantum systems whose classical limit is (in some sense) chaotic. In this subject computer 
experiments play an important role. For integrable systems the eigenvalues and eigenfunctions 
can be determined either explicitly or as solutions of simple equations. In contrast, for chaotic 
systems there are no explicit formulae for eigenvalues and eigenfunctions such that numerical 
methods have to be used. In many cases numerical observations have lead to the formulation of 
important conjectures. Such numerical computations also allow us to test analytical results which 
have been derived under certain assumptions or by using approximations. 

An important class of systems for the study of classical chaos are area-preserving maps as 
several types of different dynamical behaviour like integrable motion, mixed dynamics, ergodicity, 
mixing or Anosov systems can be found. We discuss the numerics for the corresponding quantum 
maps and illustrate some of the methods and results using the standard map and the perturbed 
cat map as prominent examples. 

Another important class of systems are classical billiards and the corresponding quantum 
billiards. In section 3 we discuss in detail the boundary integral method, which is one of the main 
methods for the solution of the Helmholtz equation, which is the time-independent Schrodinger 
equation for these systems. 



2 Area preserving maps 
2.1 Some examples 

We will restrict ourselves to area-preserving maps on the two-torus 

P ■ j2 _^ j2 ^ 
(q,p) i-> (q',p) , (2) 

where T 2 ~ R 2 /Z 2 , i.e. the map is defined on a square with opposite sides identified. The 
requirement that the map P is area-preserving is equivalent to the condition that det DP = 1, 
where DP is the linearization of the map P. The natural invariant measure on T 2 is the Lebesgue 
measure d/i = dqdp. 

As a first example let us consider the so-called standard map, defined by 

( *\ = ( q + v-f^H^q) \ modl f3 x 
V p' J V P~-t sin ( 27r <?) / 

One easily checks that this map is area-preserving. Fig. [l] shows some orbits (i.e. for different 
initial points (q,p) the points (q n ,Pn) = P n (q,p) are plotted for n < 1000) of the standard map 
for different parameters k. For k = an initial point (q,p) stays on the horizontal line and in q it 
rotates with frequency p. So for irrational p the corresponding line is filled densely. For n > 0, the 
lines with rational p break up into an island-chain structure composed of (initially) stable orbits 
and their corresponding unstable (hyperbolic) partner. For small enough perturbation there are 
invariant (Kolmogorov-Arnold-Moser or short KAM) curves which are absolute barriers to the 
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Figure 1: Examples of orbits in the standard map for different parameters k. 

motion (for a more detailed discussion of these aspects the review |TJ is a good starting point). 
For stronger perturbations, e.g. k = 1 or k = 1.5, the stochastic bands become larger and for even 
stronger perturbation (see the picture for k = 3.0) there appears to be just one quite big stochastic 
region together with the elliptic island. The elliptic islands coexist with regions of irregular motion, 
therefore the standard map is an example of a so-called system with mixed phase space or, more 
briefly, a mixed system. Whether the motion in those stochastic regions is ergodic is one of the 
big unsolved problems, see |2| for a review on the coexistence problem. For some recent results 
on the classical dynamics of the standard map, in particular at large parameters, see 0-0- 

An alternative way to specify a map P : T 2 — > T 2 is to use a generating function S(q', q), from 
which the map is obtained by 

^ = dS(q',q) = dS(q',q) 

dq dq' 
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One easily checks that 



1 K 

S(q', q) = -Aq ~ q'f + 7-5 cos(2vrg) 



(5) 



is a generating function for the standard map @. 
Another important class are perturbed cat maps 



like 
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(7) 



is a matrix with integer entries (ensuring the continuity of the map), det A = 1 (area preservation) 
and Tr A > 2 (hyperbolicity). The perturbation G(q) is a smooth periodic function on [0, 1[. For 
k = the mapping is Anosov (see e.g. ||), in particular it is ergodic and mixing. Moreover, 
following from the the Anosov theorem the map (^) is structurally stable, i.e. it stays Anosov as 
long 



y/(Tr A) 2 — 4 - Tr A + 2 
2max g |G"(g)|v / rTAi; 
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in particular the orbits are topologically conjugate to those of the unperturbed cat map. For 
larger parameters there are typically elliptic islands, so it becomes a mixed system. 
A common choice for A and the perturbation is 



(9) 
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Figure 2: Examples of orbits in the perturbed cat map (0) for k = 0.3 and k = 6.5. 



For k < K max = (y/3 — l)/\/5 = 0.33... the map is Anosov. The corresponding generating 
function is given by 

S (q, q) = q' 2 ~ qq +q 2 + £^sm(27rq) . (10) 

In fig. 0a) one orbit for 20 000 iterations for the perturbed cat map @ with k = 0.3 is shown. 
The orbit appears to fill the torus in a uniform way, as it has to be asymptotically for almost all 
initial conditions because of the ergodicity of the map. For k = 6.5 fig. |2|b) shows one orbit (20 000 
iterates) in the irregular component and some orbits (1000 iterations) in the elliptic islands. 



2.2 Quantization of area— preserving maps 

For the quantization of area-preserving maps exist several approaches, see for example |9|-[T3|. f\. 



HHTBH; a detailed account can be found in [I7|, and |18| provides a pedagogical introduction to 



the subject. First one has to find a suitable Hilbert space which incorporates the topology of the 
torus T 2 , i.e. the eigenf unctions in position and momentum have to fulfil 

^(q + j)=e^(q) ; j G N (11) 
$(p + k) = e-^ k9l $(p) ; k G N . (12) 

These conditions imply that Planck's constant h can only take the values h = whh N e N. 
Thus the semiclassical limit H — >• corresponds to iV — > oo. The phases (#i,#2) G [0, 1[ 2 are at 
first arbitrary; for 9\ = 82 = one obtains periodic boundary conditions. For each N one has a 
Hilbert space Hn of finite dimension N. Observables / G C°°(T 2 ) can be quantized analogous to 
the Weyl quantization to give an operator Op(/) on 7i n- Finally, a quantum map is a sequence 
of unitary operators Un, N e N on a Hilbert space T~Cn- The quantum map is a quantization of a 
classical map P on T 2 , if the so-called E gorov property is fulfilled, i.e. 

Jim ||^ 1 Op(/)f/ JV -Op(/oP)|| = V/GC°°(T 2 ) . (13) 

This means that semiclassically quantum time evolution and classical time evolution commute. 

So the aim is to find for a given classical map a corresponding sequence of unitary operators. 
Unfortunately, this is not as straight forward as the quantization of Hamiltonian systems and a 
lot of information on this can be found in the above cited literature and references therein. One 
of the simplest approaches to determine U ^ corresponding to a given area-preserving map uses 
its generating function to define 



{U N )j,j := {qj>\U N \qj) 



1 



d 2 S(q',q) 



dq'dq 



1/2 

exp(2mNS(q f , q] )) (14) 

q'=Qji,q=qj 



with qj = j/N, qjr = j'/N, where j,j' = 0, 1, . . . ,N — 1. In the same way one may (and for certain 
maps which cannot represented in terms of S(q', q) one has to) use other generating functions 
such as S(p',p) or S(q,p); usually these will lead to different eigenvalues and eigenf unctions. The 
question is to determine conditions on the generating function S(q r , q) such that Un is unitary 



and fulfils the Egorov property fll3|). To my knowledge this question has not yet been fully 
explored, even though the quantum maps studied in the literature provide both examples and 
counterexamples. We will leave this as an interesting open question. 
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For the examples introduced before the quantization via (II) can be used. For the standard 
map we get 



1 



exp 



m . 2 mN 



2tt 
<■<>« I ^ 



(15) 



with j,j' = 0,...,N — 1. A quantization of the standard map which takes the symmetries into 
account can be found in |19| . For the perturbed cat map (||) one gets using its generating function 



\ exp (^ (j ' 2 " j,j + j2 + K//(47r2) sin ( 27r -?/A0) 



For the unitary operator one has to solve the eigenvalue problem 

U N 1p n = X n tp n with 71 = 0, . . . , N - 1, ^ n G 



(16) 



(17) 



Here X n is the n-th eigenvalue and the corresponding eigenvector ip n consists of N complex compo- 
nents, where N is the size of the unitary matrix Un- Because of the unitarity of Un the eigenvalues 
lie on the unit circle, i.e. |A n | = 1. 

Let us discuss some of the numerical aspects relevant for finding the solutions of (|T?D without 



going into implementation specific details (see the appendix and [2(| for an implementation using 
Python). 

Computing the eigenvalues of ([17]) consists of two main steps 

• Setting up the matrix Un- 

The computational effort increases proportional to N 2 (unless each matrix element requires 
further loops) as we have to fill the N 2 matrix elements. 

The memory requirement to store Un is 16 N 2 Bytes (for a IEEE-compliant machine a double 
precision floating point number requires 8 Bytes; as we have both real and imaginary part 
we end up with 16 Bytes per matrix element). 

• Computing the eigenvalues: 

The computational effort for the matrix diagonalization (typically) scales like N 3 . 
Usually one will use a black-box routine such as one from the NAG-library |BT| or from 



LAPACK |22[ . To my knowledge there are no routines which make use of the fact that the 
matrix Un is unitary so we may for example use the NAG routine F02GBF or the LAPACK 
routine ZGEES (or the more recent routine ZGEEV which is faster for larger matrices, e.g. 
iV > 500) which compute all eigenvalues of a complex matrix. 

For certain maps specific optimizations are possible, see e.g. [W\ for the standard map. For 
this type of mapping a different approach employing a combination of fast Fourier transform 
and Lanczos method reduces the computational effort to iV 2 lniV [23|. 

After successful compilation and running of the program it is useful to see whether the eigen- 
values really lie on the unit circle. In fig. |3] this is illustrated for N = 200 and the standard map 
with k = 1.5. For small iV the running times of the program for setting up the matrix Un and its 
diagonalization is just a matter of minutes. For example on an Intel Pentium III processor with 
666 MHz one needs just 6 minutes to compute the eigenvalues of ( |TT| ) when N = 1000. However, 
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Figure 3: Plotting the eigenvalues of Un allows to check the numerical implemen- 
tation and the unitarity of Un', the picture shows for iV = 200 and k = 1.5 the 
eigenvalues A n for the quantized standard map (|15|). 



for N = 3000 already 140 MB of RAM are required to store Un and the computing time increases 
to 6 hours. Depending on available memory, computing power, patience and motivation one may 
use larger values of N. 

Let us conclude this part with a more technical remark: In addition to the choice of com- 
puter language, compiler, optimizations and algorithm there is one very important component 
for achieving good performance when doing numerical linear algebra computations: the BLAS 
(Basic Linear Algebra Subprograms). Libraries such as LAPACK defer all the basics tasks like 
adding vectors, vector-matrix or matrix-matrix multiplication to the BLAS such that highly op- 
timized (machine-specific) BLAS routines should be used. Most hardware vendors provide these 
(of differing quality). Recently the software system ATLAS (Automatically Tuned Linear Algebra 
Software) [E3] has been introduced which generates a machine dependent optimized BLAS library. 
For some computers ATLAS-based BLAS can be even faster than the vendor supplied ones! 



2.3 Eigenvalue statistics 

One central research line in quantum chaos is the investigation of spectral statistics. It has been 
conjectured p5[ that for generic chaotic systems the eigenvalue statistics can be described by 
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Figure 4: (a) Level spacing distribution P(s) and (b) cumulative level spacing dis- 
tribution J(s) for the perturbed cat map (Q) with k = 0.3 and iV = 3001. 
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Figure 5: (a) Level spacing distribution P(s) and (b) cumulative level spacing dis- 
tribution I(s) for the perturbed cat map (p2|) with = 0.012 and N = 3001. 
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random matrix theory, whereas generic integrable systems should follow Poissonian statistics ||26|| . 
To study the eigenvalue statistics for quantum maps one considers the eigenphases (p n G [0,27r[, 
defined by A n = e llfn (in the following we will also call (p n levels in analogy to the energy levels 
for the Schrodinger equation). The simplest statistics is the nearest neighbour level spacing 
distribution P(s) which is the distribution of the spacings 

N 

s n ■= TT^ri+i ~ <p n ) with 71 = 0, . . . , JV - 1 and tp N := y? . 

271 

The factor j- ensures that the average of all spacings s n is 1. To compute the distribution 
practically one chooses a division of the interval [0, 10] (usually this interval is sufficient, but more 
precisely the upper limit is determined by the largest s n ) into b bins and determines the fraction 
of spacings s n falling into the corresponding bins. If iV is too small it is better to consider instead 
of P(s) the corresponding cumulative distribution 

/« := #< " (18) 

which avoids the binning and results in a smoother curve. 

Fig. H shows for the perturbed cat map @ with k = 0.3 the level spacing distribution P(s) and 
the cumulative level spacing distribution I(s) for N = 3001. For this parameter value k the map is 
still Anosov so one expects that the correlations of the eigenphases follow random matrix theory; 
in particular because the perturbation should break up the number theoretical degeneracies which 
lead to non-generic spectral statistics for the cat maps at k = [p7,^8fl. In [^, it is shown 



that for all perturbations which are just a shear in position one of the symmetries of the cat map 
survives, so that the statistics are expected to be described by the circular orthogonal ensemble 
(COE). In the limit N — > oo this is the same as the Gaussian orthogonal ensemble (GOE). In 
fig. £| we show the Wigner distribution Pwigner(s) which is very close to the COE distribution, 



7T ( IX 2 



Pcue(s)^— s 2 exp(-- S 2 ) (20) 



^COE(s) W iwigner(s) = -Sexp (--S ) . (19) 

and for comparison the CUE distribution 

32 2 / 4 
— s exp 

7T Z \ 7T 

and the Poisson distribution (expected for generic integrable systems) 

J Ppoisson(s) = e _s . (21) 

The agreement with the expected COE distribution is very good. 

A specific example, which breaks the above mentioned unitary symmetry and thus leads to 
CUE statistics, uses two shears, one in position and one in momentum 



where 



q p ,)=(A°P,oP p )( q p ) (22) 



1 = I 41 24 ) {2:; ' 



and Pq(q,p) = (q+ K q G(p),p), P p (q,p) = (q, p + K p F(q)) with F(q) = ±(sm(2nq) - sin(47rg)) and 
G{p) = 2^:(sin(47rg) — sin(27rg)). For the corresponding quantum map with k p = K q = 0.012 and 
iV = 3001 the level spacing distribution is shown in fig. |||. One observes very good agreement 
with the CUE distribution. 
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2.4 Eigenfunctions 



Another interesting question concerns the statistical behaviour of eigenfunctions, and more specifi- 
cally for quantum maps the eigenvector statistics and the properties of phase space representations 
like the Husimi function. 



2.4.1 Eigenvector distributions 

Consider an eigenvector ip of a quantum map given by the N numbers Cj G C, j = 0, N — 1. 
The distribution P{ip) is given (similarly to the level spacing distribution) by 



b 

^#{a < |c,f < b} = j P(^)diP 



(24) 



Let us first discuss the corresponding random matrix results (see e.g. [31, 32]). For the COE the 
eigenvectors can be chosen to be real and the coefficients Cj, j = 0, . . . , N — 1, only have to obey 
the normalization condition 

N-l 

^2^ = 1 with cj G R . (25) 

3=0 

Thus the joint probability for an eigenvector c = (cq, . . . , cn-i) G M. n is 



where the prefactor ensures normalization. So the probability of one component to have a specific 
value y is given by integrating P^ OE (c) over all other components, 



P% OE (y) = j s(y - cl)p^(c) dc ■ ■ ■ dc N ^ = -L T *Wf )/2) (i - y) iN ~ 3)/2 ■ 

The mean of P^ OE (y) is L yP^ OE (y) dy = 1/N. So using the rescaling rj = yN gives 



(27) 



^ 0E <"> = 7^,fW%) 11 - * ,Nr ~ m ' (28) 



In the limit of large N one gets the so-called Porter-Thomas distribution |33| 

^ OE (^) = ^exp(-r / /2) , (29) 

y 

and the corresponding cumulative distribution I(y) — f P(y') dy' reads 

o 

I{r,) = erf ( yfijt) . (30) 
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Figure 6: (a) Eigenvector distribution for the perturbed cat map (P) with iV = 1597 
and k = 0.3. In comparison with the asymptotic COE distribution, eq. (p8|), dashed 
line. The inset shows the same curves in a log-normal plot. In (b) the cumulative 
distribution is shown and in the inset a plot of the absolute value of the components 
Cj, j = 0, . . . , iV — 1 of the corresponding eigenvector ip n= 2o is displayed. 
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Figure 7: (a) Eigenvector distribution of an eigenvector for the perturbed cat map 
(E3) with iV = 1597 and = 0.012 is shown in comparison with the asymptotic 

CUE distribution, eq. (^), dashed line. The inset shows the same curves in a log- 
normal plot. In (b) the corresponding cumulative distributions are shown and in 

In) 

the inset a plot of the absolute value of the components Ca , j = 0, . . . , N — 1 of the 
corresponding eigenvector ipn=2 is displayed. 
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Fig. ^| shows an example for the eigenvector distribution of an eigenstate of the perturbed cat map 
(H) with k = 0.3 and N = 1597. There is good agreement with the expected COE distribution, 
eq. (H), shown as dashed line. 



Finally, let us consider again the map (22) which shows CUE level statistics. From this one 
would expect that also the eigenvector statistics follows the CUE. Similar to the case of the COE 
one has the normalization condition 



JV-l 
3=0 



with c~ G C 



(31) 



Thus the joint probability for an eigenvector c = (co, . . . , cn-i) G C reads 

N-l 
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The probability of one component to have a specific value y is given by integrating P 1 
all other (complex) components, 
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Again as for the COE, the mean of P^j VE (y) is 1/N and the rescaling t] := yN leads to 
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Figure 8: Eigenvector distribution of an eigenvector for the unperturbed (i.e. k — 0) 
cat map (|9|) with iV = 1597. This is compared with the asymptotic semicircle 
law, eq. (|37|) . The inset shows the corresponding eigenvector (compare with the 
eigenvectors shown in the previous two figures). 
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which has mean 1. In the large X limit we have 



-.CUE 



(77) = exp(-77) 



(35) 



and the cumulative distribution simply is 

I CVE (rj) = 1 - exp(-77) 



(36) 



Fig. shows P{rf) for one eigenvector of the perturbed cat map fl22[). There is good agreement 
with P CUE (r/). 

A different distribution is obtained for unperturbed cat maps: for certain subsequences of 
prime numbers (which depend on the map) the distribution of 77 = |Re?/> tends to the semicircle 
law, 



for 7] < 1 
for 7] > 1 



(37) 



sec 



34| for details (see also ||35|| ). In fig. |8| we show an example of an eigenstate with A" = 1597 



for the quantum map corresponding to the map (^) with k = 0. For this X the map fulfils the 
conditions of MM and one observes a nice semicircle distribution of the eigenvector. However, it 
seems that the approach to the asymptotic distribution is slower than for the case of the random 
matrix situations. 



2.4.2 Husimi functions 

A different representation of eigenstates is to consider a phase space representation, like for ex- 
ample the Husimi function, which allows for a more direct comparison with the structures for the 
classical map. Without going into the mathematical details, the Husimi representation is obtained 
by projecting the eigenstate onto a coherent state centred in a point (q,p) G T 2 , 



H n (q,p) = \(C q M n )\ 



N-1 



N-1 



^{C q , P \qj){qj 



j=0 



N-1 



y^(^g,p|gj) c j 



3=0 



^(2A0 1/4 exp (-7cN{q 2 - ipq)) exp{7cX{-q] + 2{q - ij%)) 
3=0 

2 



t)a I it: X [qj-^-q + ip 



iX ) Cj 



Here qj = ^(^2 + j), j — 0, . . . ,N — 1 and ^(Zjr) is the Jacobi-Theta function, 



HZ\r) = E 



iirTn 2 +2mZ 



with Z, t e C, Im(r) > . 



(3? 



(39) 



The coefficients Cj are the components of the eigenvector ip n in the position representation as 
obtained from the diagonalization of Un (for other generating functions than the one used in 
eq. ( [Lip one has to adapt eq. (p8|)). 

If one wants to compute a Husimi function on a grid of X x X points on T 2 the computational 
effort grows with A^ 3 . So for computing all Husimi function of a quantum map for a given A^ the 
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Figure 9: In a) a Husimi function H n (q,p) of the perturbed cat map @ with « = 0.3 
is plotted which shows the expected 'uniform' distribution. Here black corresponds 
to large values of H n (q,p). In b) for k = 6.5 a state localizing on one of the elliptic 
islands is shown (compare with fig. |||). 
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Figure 10: Examples of Husimi functions for the standard map with k = 3.0 and 
N = 1600. 
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Figure 11: Examples of Husimi functions H n (q,p) for the standard map with k = 1.5 
and iV = 1600 for n = 0, . . . , 19. (Compare with fig. []].) 
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computational effort grows with iV 4 . Already for moderate N this can be quite time-consuming, 
but even more importantly, usually one also wants to store all these Husimi functions on the hard- 
disk which limits the accessible range of N. Sometimes a smaller grid, e.g. of size lOv^A x lOy/N 
can be sufficient which reduces the growth of the computational effort to iV 2 for a single Husimi 
function and to N 3 for all Husimi functions at a given N. Even then one still needs 800 N 2 Bytes 
to store these on the hard-disk. For example for N = 1600 this roughly leads to 2 GB of data and 
for iV = 3000 one needs approximately 7 GB. However, there are also cases where a finer grid, 
e.g. 2N x 2N is necessary. 

Theoretically one expects that for iV — > oo the Husimi functions concentrate on those regions 
in phase space which are invariant under the map (this follows from the Egorov property). So for 
ergodic systems the expectation is that (in the weak sense) 

H n {c[,p) —> 1 with n — 0, N — 1 as the matrix size N — ► oo . (40) 

The precise formulation of this statement is the contents of the quantum ergodicity theorem for 



maps |36| (see [3^J for the case of discontinuous maps). The quantum ergodicity theorem only 
makes a statement about a subsequence of density one (i.e. almost all states) which for example 
leaves space for scars, i.e. eigenstates localized on unstable periodic orbits. For systems with mixed 
phase space one (asymptotically) expects localization in the stochastic region(s) and on the tori 
in the elliptic regions. 

In fig. |^a) we show for the perturbed cat map with k = 0.3 the Husimi function for the same 
eigenstate as in fig. || As expected it shows a quite uniform distribution (of course with the usual 
fluctuations). In contrast for k = 6.5 there are eigenstates such as the one shown in fig. ^jp) which 
localizes on the elliptic island (compare with fig. 0). 

In some sense more interesting are the Husimi functions for mixed systems as the classical 



dynamics shows more structure. In fig. [10| we show some examples for the standard map with 
k = 3.0 Fig. |T0"|a) shows a Husimi function which is spread out in the irregular component. In 
contrast in b) the Husimi function localizes on a torus around the elliptic fixpoint. The Husimi 
function in c) shows quite strong localization around the small elliptic island of a periodic orbit 
with period 4. This island is so small that it is not visible in fig. |l|. Therefore, the Husimi 
function displayed in fig. |TD|d) indicates that the region of 'influence' of this island is much larger 
than the area of the island. This region is also visible in the Husimi function in fig. |10|a), as the 
irregular state has a very small probability in the regions around these islands. A longer sequence 



of Husimi functions for the standard map with k = 1.5 shown in fig. |Tl] illustrates the different 
types of localized states (compare with fig. [I]) . 



3 Billiards 

3.1 Classical billiards 

A two-dimensional Euclidean billiard is given by the free motion of a point particle in some 
domain Q C M 2 with elastic reflections at the boundary dQ. Depending on the boundary one 
obtains completely different dynamical behaviour, see fig. [12] where this is illustrated by showing 
orbits of billiards in a circle, a square and an ellipse, which are all integrable giving rise to regular 
motion. In contrast the Sinai billiard (motion in a square with a circular scatterer), the stadium 
billiard (two semicircles joined by parallel straight lines) and the cardioid billiard show strongly 
chaotic motion (they are all proven to be hyperbolic, ergodic, mixing and if-systems) . 
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Integrable systems 




Figure 12: Billiard dynamics in integrable and chaotic systems. 

As the motion inside the billiard is on straight lines it is convenient to use the boundary to 
define a Poincare section, 

V:={(s,p)\se[0,\dn\],pe[-l,l]} . (41) 

Here s is the arclength along <9f2 and p = (v,T(s)) is the projection of the unit velocity vector 
v after the reflection on the unit tangent vector T(s) in the point s e d£l. The Poincare map is 
then given by 

P : V -> V 

i.e. for a given point £ = (s,p) one considers the ray starting at the point r(s) e <9f2 in the direction 
specified by p and determines the first intersection with the boundary, leading to the new point 
£' = (s',p'). Explicitly, the Cartesian components of the unit velocity v of a point particle starting 
on d£l at r(s) are determined by the angle (3 G [— 7r/2,7r/2] measured with respect to the inward 
pointing normal N = (—T y ,T x ). The velocity in the T,N coordinate system is denoted by 
(p,n) = (sin f3, cos f3), so that in Cartesian coordinates 

v = (v x , v y ) =^ ^ (p, n) = (r x p + N x ^/l^, T y p + N y y/T^j . (42) 
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Starting in the point r(s) G dQ in the direction v, the ray r + tv intersects dQ at some point 
r' = (x',y'). If the boundary is determined by the implicit equation 

F(x,y) = , (43) 

the new point r' can be determined by solving 

F(x + tv x ,y + tv v ) = . (44) 

For non-convex billiards there are points ^ = (s,j)) £ ? for which there is more than one solution 
(apart from t = 0); obviously the one with the smallest t > has to be chosen. The condition 
can be used to remove the t = solution analytically from (H5). If F is a polynomial in x 



and y this allows to reduce the order of fl44 ) by one. This approach has for example been used for 
the cardioid billiard leading to a cubic equation for t, see |3S[] for details. From the solution t one 
gets the coordinates (x 1 , y') = (x, y) + tv which have to be converted (in a system dependent way) 
to the arclength coordinate s' (in many practical applications there is a more suitable internal 
variable, for example the polar angle etc.). The corresponding new projection of the momentum 
is given by p' = —(v,T(s')). 

3.2 Quantum billiards 

For a classical billiard system the associated quantum billiard is given by the stationary Schrodinger 
equation (in units h = 2m = 1) 

-Aip n {q) = E n ip n (q) , qett (45) 

with (for example) Dirichlet boundary conditions, i.e. if) n (q) = for q G dQ. Here A denotes the 
Laplace operator, which reads in two dimensions 



\dqf dq 2 



A =\ ^ + ^2) • ( 46 ) 



In the Schrodinger representation the state of a particle is described in configuration space by 
a wave function ip G L 2 (Q), where L 2 (Q) is the Hilbert space of square integrable functions on Q. 
The interpretation of ip is that J D \ip(q)\ 2 d 2 q is the probability of finding the particle inside the 
domain D C Q. 

Due to the compactness of Q, the quantal energy spectrum {E n } is purely discrete and can be 
ordered as < E\ < E2 < E% < . . . . The eigenfunctions can be chosen to be real and to form an 
orthonormal basis of L 2 (Q), 

(ip n \i>m) ■= / i>n{q)i>m{q) d 2 g = s mn . 



The mathematical problem defined by eq. fl4"5] ) is the well-known eigenvalue problem of the 
Helmholtz equation, which for example also describes a vibrating membrane or flat microwave 
cavities. For some simple domains Q it is possible to solve eq. ( |4"5] ) analytically. For example for 
the billiard in a rectangle with sides a and b the (non-normalized) eigenfunctions are given by 
i^m^iq) = sin (7mi gi /a) sin(7rn2<2 , 2/&) with corresponding eigenvalues E nijTl2 = n 2 {n\/a 2 + n\/b 2 ) 
and (711,112) G N 2 . For the billiard in a circle the eigenfunctions are given in polar-coordinates 
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Figure 13: Plot of A^ uc (x) for the stadium billiard (a = 1.8, odd-odd symmetry) to- 
gether with the contribution from the bouncing ball orbits, dashed line, see eq. fl52|). 
In b) the fluctuating part after subtraction of the contribution of the bouncing ball 
orbits is shown. 
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Figure 14: Detection of missing levels using the 5 n -statistics. Shown is Na uc (x) — 
^osc( x ) f° r the case where one level has been removed from the spectrum. One 
clearly sees the jump in the fluctuations which allows to roughly locate the place 
where the level is missing. 



by ip mn (r,{p) = J m (jmn r ) exp(imy?), where j mn is the n-th zero of the Bessel function J m (x) and 
m G Z, n G N. However, in general no analytical solutions of eq. (551) exist so that numerical 
methods have to be used to compute eigenvalues and eigenfunctions. 
The spectral staircase function N(E) (integrated level density) 



N(E) := #{n | E n < E} 



(47) 



counts the number of energy levels E n below a given energy E. N(E) can be separated into a 
mean smooth part N(E) and a fluctuating part 



N(E) = N(E) + N {luc (E) 



For two-dimensional billiards, N(E) is given by the generalized Weyl formula |39 



N(E) 



A 

47T J 



E 



47T 



E + C + ... 



(48) 



(49) 



where A denotes the area of the billiard, and C := L~ — C + , where Cr and C + are the lengths 
of the pieces of the boundary dQ with Dirichlet and Neumann boundary conditions, respectively. 
The constant C takes curvature and corner corrections into account. 

The simplest quantity is the <5 n -statistics, which is obtained from the fluctuating part of the 
spectral staircase evaluated at the unfolded energy eigenvalues x n := N(E n ) 

S n := N (E n ) - N(E n ) = n -\-x n , (50) 
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where 



The quantity 5 n is a good measure for the completeness of a given energy spectrum. For a complete 
spectrum S n , or equivalently N auc (x), should fluctuate around zero. Fig. |T3|a) shows Nq uc (x) for the 
stadium billiard, which indeed fluctuates around zero. In addition there is an overall modulation 
of iVfl uc (x) which is caused by the bouncing ball orbits. They form a one parameter family of 
periodic orbits having perpendicular reflections at the two parallel walls (of length a, see fig. [l5|) 
of the stadium. The contribution of these orbits to the spectral staircase function reads [EH] 



= I E V^E? 8 (VI - v^) - (±B - ±Ve) 

71=1 ^ ' 

a i *^ 1 / / — 3n\ 
= —= > — cos 2anVE - — , 

2v K i «2 V 4 / 

v 71=1 N ' 



(52) 
(53) 



where E^ h = n 2 n 2 are the eigenvalues of a particle in a one- dimensional box of length 1, and is 
the Heaviside step function. Subtracting N^ c (x) from A^ uc (x) removes the additional oscillation, 
see fig. p~3[b) . If an eigenvalue is missing this is clearly visible by a 'jump' of S n in comparison to 



points fluctuating around 0, see fig. [14] for an example where one eigenvalue has been removed 
'by-hand'. Clearly, the energy interval in which a level is missing can be estimated from the plot. 

In the same way as for quantum maps one can study the level spacing distribution and more 
complicated statistics, like the number variance, n-point correlation functions etc., see for example 
iT|.|42]| for some further examples for the cardioid billiard. 



3.3 Computing eigenvalues and eigenfunctions for quantum billiards 

There exist several numerical methods to solve the Helmholtz equation inside a domain Q C M 2 , 

AiP(q) + k 2 i)(q) = , q G Vl\dn , (54) 
with Dirichlet boundary conditions 

V>(g) = , qedtt. (55) 



see |j43[ , which however does 
Additionally, in the context of 
for a detailed description of the 



For a good review on the determination of the eigenvalues of 
not cover finite element methods or boundary integral methods 
quantum chaos the plane wave decomposition Q| (see also [45 
method), the scattering approach, see e.g. |46|-48[[, and more recently the scaling method |49|, are 
commonly used. 

Here I will give a sketch of the derivation of the boundary integral method and discuss in 
more detail the numerical implementation. The boundary integral method reduces the problem 



of solving the two-dimensional Helmholtz equation (54) to a one-dimensional integral equation 



see e.g. 



mm 



and references therein. Of course, the general approach also applies to higher 
dimensions but we will only discuss the two-dimensional case. For studies of three-dimensional 
systems by various methods see e.g. [po}-[7D|]. Boundary integral methods are also used in many 
other areas so that it is impossible to give a full account. For example they are also commonly used 
in acoustics, see e.g. [|n]] and the detailed list of references therein. Finally, the boundary integral 
method provides a starting point to derive the Gutzwiller trace formula, see e.g. W% 
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3.3.1 Boundary integral equation 

Let G(q, q') be a Green function of the inhomogeneous equation, i.e. 

(A + k 2 )G k (q,q')=5(q-q') . 
Considering the integral over Q of the difference ip(q')-(ffi§—G k (q, one obtains 

J mq')A'G k (q,q') -G k (q,q')A'i;(q')} d 2 q' = J ^)5{q - q') d 2 q> . 



(56) 



(57) 



n n 
Using the second Green theorem gives the Helmholtz representation 



dG 



M)-^(q,q')-G k (q,q')^(q' 



dn' 



ds' 







q e Q \ dtt 

qedn 

else 



(58) 



Here q' = q{s') and = n(s')V with n(s) = (q' 2 (s) , —q[(s)) denoting the outward pointing 
normal vector, where (qi(s), ^(-s)) is a parametrization of the billiard boundary dfl in terms of 
the arclength s, oriented counterclockwise. Special care has to be taken to obtain the result for 
q G dQ, see e.g. []5T],[F4]]. (When q is in a corner of the billiard the factor | has to be replaced by 
where 9 is the (inner) angle of the corner.) For Dirichlet boundary conditions one obtains 



j u{s')G k {q,q') ds' = , q E dtt , 



(59) 



where 



<)<1 







u(s) := —ip(q(s)) := n(s)Vi/)(q(s)) := n(s) hm Vipiq ) 
on q'^q(s) 

q'en\dn 



(60) 



is the normal derivative function of i/j. 

In two dimensions a Green function for a free particle is given by the Hankel function of first 
kind 



G k (q,q') 



1 - r W 1 



- (k \q -q'\) = -- [J (k \q-q'\)+i Y (k \q - q'\ 



(61) 



Since H^\z) ~ - In z for z — > 0, the Green function G k (q, q') diverges logarithmically such that 
it is more convenient to derive an integral equation whose kernel is free of this singularity. To 
that end one (formally) applies the normal derivative to eq. (|58|). More carefully one has to 



consider a jump relation for the normal derivative function, see e.g. |51"|, [T4|| . The result is 

u{s) = -lj ^-G k (q(s), q(s')) u(s') ds' . (62) 



an 

For the derivative of the Green function one obtains 



^-G k (q(s), q(s')) = ^ cos(0( S , s')) (k r(s, s')) 



(63) 
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where r(s, s') = |qr(s) — q(s')\ is the Euclidean distance between the two points on the boundary 
and 

C0S ^ M= rMm^m . (64) 

r(s, s') 

This gives the integral equation for the normal derivative u(s) 

u(s)= J>Q k (s,s')u(s')ds' , (65) 
an 

with integral kernel 

Q k (s,s') = -^cos0( S , S ') (kr{s,s')) . (66) 
Eq. (|65D is a Fredholm equation of second kind which has non-trivial solutions if the determinant 

D(k) := det(l - Q h ) (67) 
vanishes. Here Q k is the integral operator on dVL defined by 

Qk{u{s)) = lQ k {s,s') u(s') ds' . (68) 



an 

For eigenvalues E n of the Helmholtz equation with Dirichlet boundary conditions one has 



D(k) = for k = \^E^ } see e.g. [0] for a detailed proof. However, for ImA; < there can be 
further zeros of D(k) which (for the interior Dirichlet problem) correspond to the outside scattering 
problem with Neumann boundary conditions ||76|-[78|1 (see also |51|). Explicitly this can be seen 



from the factorization D(k) = D(0)D int (k)D cxt (k), where the factors can be written exclusively 
in terms of the interior and exterior problem. More aspects concerning the additional spurious 
solutions will be discussed in sec. |3.3.5[ 



Before turning to the numerical implementation, let us discuss the behaviour of the integral 
kernel for small arguments. The Hankel function H± (x) reads for small arguments 

h[ 1] (k t(s, s')) rr^ r , for s - s' -> . (69) 

nk\s — s \ 

This singularity is compensated by the behaviour of 

cos0(s, s') ~ — -/t(s) |s — s'| , for s' — > s , (70) 

where k(s) is the curvature of the boundary in the point s. Here the curvature is defined by 
K ( s ) = Qii^I'li. 3 ) ~ <?2( s )<?i( s ) sucn ^at for example k(s) = 1 for a circle of radius one. Thus for 
the integral kernel we obtain 

Q k (s, s') -> ^-k(s) , for s - s' -> . (71) 

Z7T 
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Figure 15: Geometry of the desymmetrized stadium billiard. 



3.3.2 Desymmetrization 

For systems with symmetries the numerical effort can be reduced by considering instead of the 



full system the symmetry reduced system with the corresponding Green function, see e.g. |56 
For a reflection symmetry with respect to the gi-axis one has 

Gt(q,q'):=G k (\q-q'\)±G k (\q-(q' 1 ,-q^\) , (72) 

where + applies to the case of even eigenfunctions (i.e. Neumann boundary conditions on the 
symmetry axis) and — to odd eigenfunctions (i.e. Dirichlet boundary conditions on the symmetry 
axis) . 

For a two-fold reflection symmetry (as in the case of the stadium billiard, see fig. To* for a 



sketch of the geometry and notations) one has altogether four different subspectra, corresponding 
to DD, DN, ND and DD boundary conditions on the symmetry axes q\ and q2, respectively. For 
example for Dirichlet-Dirichlet boundary conditions on the q\- and flaxes the Green function 
reads 

GT = G k (\q -q'\)- G k (\q - (q[, -q' 2 )\) + G k (\q - (-q[, -q' 2 )\) - G k (\q - (-q[, q' 2 )\) . (73) 

For Neuman boundary conditions on these two axes one gets 

GT = G k (\q -q'\) + G k (\q - (q[, -q' 2 )\) + G k (\q - (-q[, -q' 2 )\) + G k (\q - (-q[, q' 2 )\) . (74) 

The advantage of exploiting the symmetries of the system is two-fold: firstly, we can separate 
the eigenvalues and eigenfunctions for the different symmetry classes, which is necessary for the 
investigation of the spectral statistics. Secondly, the numerical effort is reduced, since the integral 
over the whole boundary dVL is reduced to an integral over a part of the boundary, which in the 
above examples is half or a quarter of the original boundary. The boundary along the symmetry 
axes need not be discretized as the boundary condition is already fulfilled by construction. Of 
course, for other geometries different choices for G can be more appropriate. 
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Figure 16: In a) the three smallest singular values are shown as a function of the 
energy E = k 2 for the stadium billiard with a = 1.8 and odd-odd symmetry. The 
eigenvalues are located at the minima of the first singular value. The second and 
third singular values allow to locate places with near degeneracies as next to k 2 = 90, 
which can be resolved by magnification of the corresponding region, see fig. pj]. In 
b) | det(Afc)| is shown. The minima tend to be not as pronounced as those of the 
singular values. 
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Figure 17: A magnification of fig. [16] shows that the singular value decomposition 
method easily allows to locate nearly degenerate energy levels. 




Figure 18: Plot of real and imaginary part of det(v4fc) as a function of fc; the eval- 
uation was done for 10 times as many points in k 2 than for fig. [T7]. Approximately 
simultaneous zeros correspond to minima of | det(^4fe)|. The locations of the eigen- 
values are marked by squares. 
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3.3.3 Finding the eigenvalues 

In the numerical computations the integral over the boundary is replaced by a Riemann sum. 
(There also exist more refined methods using polynomial approximations combined with GauB- 
Legendre integration, see e.g. f60f , which allow for a less fine discretization.) Let As = C/N be 



the discretization length of the boundary of length C into N pieces. Then we have 

N-l 

u(sj) = As Qkjsj, Sj) u(sj) , (75) 

where Sj = (i + 1/2) As, % = 0, . . . , N — 1. Equation (|75|) can be written in matrix form as 

A k u = , with Aij = 5ij - As Q k {si, sj) . (76) 

Recall that for Sj = Sj the kernel Qk(si,Sj) reduces to the result given in eq. (|7l|). The solutions 
of this linear equation provide approximations to the eigenvalues k\ and eigenvectors u n . This 
leads to the problem of finding the real zeroes of the determinant 

det(A fc ) = (77) 

as a function of k = y/E, where A k is a dense, complex non-Hermitean matrix. Due to the 
discretization of the integral the determinant det(A^) will not become zero but only close to zero 
(actually, the discretization shifts the zeros slightly away from the real axis, see ]58"|, CT] ) . 

In the numerical computations it is very useful ]60| to compute the singular values of the 
matrix A instead of its determinant. The singular value decomposition of a complex matrix is 
given by the product of an unitary matrix U, a diagonal matrix S and a second unitary matrix V 

A = USV^ . (78) 

The diagonal matrix S contains as entries SV« the singular values of A and we have | det A\ = 
]~|SVj|. Since the original integral equation has been discretized, the smallest singular value in 
general never gets zero, but just very small, see fig. |TB|. Thus the minima of the smallest singular 
value provide approximations to the eigenvalues of the integral equation. For the numerical 
computation of the singular value decomposition one may for example use the NAG routine 
F02XEF or the LAPACK routines ZGESVD or ZGESDD. It turns out that the (more recent) routine 
ZGESDD is significantly faster (factor 3-5, at the expense of a higher memory consumption), in 
particular when also singular vectors are computed. 

The advantage of the singular value decomposition in comparison to locating the zeros of the 
determinant is that degeneracies of eigenvalues can be detected by looking at the second singular 
value, which also gets small when there are two nearby eigenvalues (similarly higher degeneracies 
can be found by looking at the next singular values). In fig. |T6|a) an example of the behaviour 
of the three smallest singular values is shown in the case of the stadium billiard (a = 1.8) with 
Dirichlet boundary conditions. For comparison a plot of |det(Afc)| is shown in fig. |T^b). One 
clearly sees that the singular value decomposition provides more information. For example, next 
to k 2 = 90 the minimum of | det(A&)| looks slightly broader than the others, however, this does 
not give a clear indication that there might be more than one eigenvalue. In contrast, the singular 



value decomposition method allows to resolve such kind of near-degeneracies efficiently, see fig. 17 



Of course, this information is also available via det(A k ), see fig. [18] where its real and imaginary 
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part are plotted separately. Here (approximately) simultaneous zeros correspond to minima of 
|det(Afc)|. However, notice that compared to the singular value decomposition approach much 
more discretization points in E = k 2 are necessary. 

To determine all energy levels in a given energy interval [Ex, E 2 ] one proceeds in the following 
way: first one computes the singular values at equidistantly chosen points k 2 G [Ex, E 2 ]; the energy 
is chosen as variable because for two-dimensional billiards the mean distance between two energy 
levels is approximately constant and according to the generalized Weyl formula ([K]) given by 
The finer the step size is chosen the easier the minima can be resolved, however, at the same time 
the computing time to cover a given energy range increases correspondingly. The actual step size 
is a compromise between these two aspects; good results have been achieved by using a step size of 
the order of | (for systems with many near level degeneracies, e.g. integrable or near-integrable 
systems, a smaller step size can be helpful). 

The matrix size N is chosen according to N = = b^-, such that one obtains b discretization 
points per units of the inverse of the de Broglie wave length A = ^ along the boundary C. Typical 
choices for b are between 5 and 12 depending on the system and the wanted accuracy. 

From the first scan one locates all minima of the smallest singular value. If also the second 
singular value has a minimum next to a minimum of the first one, one has to use a refined 
discretization in E around the minimum (the numerical implementation is a bit more sophisticated, 
in order to account for several special situations, so that only a minimal number of additional points 
need to be computed). Once an isolated minimum is found, an approximation to the eigenvalue 
can be computed by different methods. Either one can perform a refined computation around 
the minimum, which can be quite time-consuming, or one can use a local approximation by a 
parabola |79|| . A linear interpolation also gives good results: From the three points l:(kf, SVi(k%)), 
2:(&f, SVxik 2 )), 3:(/c|, SVx{k 2 )), characterizing a minimum of the first singular value, one has two 




Figure 19: Magnification of fig. [IT] around the minimum with k 2 m 96.5 for different 
matrix sizes N. One nicely sees the pronounced parabolic structure for N = 45 
which gets smaller for larger N. 
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different lines 12 and 23 with different slopes, of which the line with the larger slope has to 
be chosen. The intersection of this line with the zero axis gives a good approximation to the 
eigenvalue, which one can refine if necessary. By repeatedly applying this for all minima, all 
energy levels in a given interval can be found. In fact, it is possible to develop a computer 
program which takes care of all this such that all levels can be found automatically. 

A good check of the completeness is provided by considering the 5 n statistics, see the example 
in sec. |3.2| . The accuracy of the computed eigenvalues can be estimated from the bracket of the 
minimum given by the three points 1,2,3 if the matrix dimension N is large enough. For N too 
small (for a given resolution in E) one does not obtain a peaked, but a broad minimum. This is 
illustrated in fig. |1£] by magnifying fig. [IT] around the minimum with k 2 ~ 96.5 for different N. 
One clearly sees the parabolic structure around the minimum for smaller N and for larger N one 
recovers the essentially linear behaviour of the smallest singular value. 

Tests of the accuracy of the method can be obtained by considering a system where the 
eigenvalues are known. For example for the circular billiard the eigenvalues can be computed with 
arbitrary accuracy. Also billiards where the eigenvalues can be computed by other methods (e.g. 
conformal mapping method 51]) allow a determination of the accuracy of the method. For 
a study of the scaling of the error for various billiards see In addition computations of the 
normal derivative function u n (s) and the eigenfunction (both inside and outside of Q) allow to 
check the quality of the numerical method and program. 



3.3.4 Computing eigenfunctions 

From a minimum of the smallest singular value we obtain an approximation of the eigenvalue 
and at the same time the corresponding singular vector u gives an approximation to the normal 
derivative function u(s). The NAG routine F02XEF scales the singular vector such that its first 
component is real. Thus for a correct solution also the other components should be essentially 
real, which provides another check for the implementation of the method and the accuracy of the 
eigenvalues. 

The eigenfunction in the interior of the domain Q can now be calculated from the normal 
derivative function, 

^(q) = - l -jH^(k\q-q(s)\)u{s)ds, for q G Q\dQ . (79) 
an 

The computation of the eigenfunction can be simplified by taking into account that 

j J (k \q - q(s)\) u{s) ds = , (80) 
an 

because the Jo-part of Gk(q,q') is a solution of the homogeneous equation corresponding to 
eq. ([SB]). Thus (|7|) is equivalent to 

^(q) = ^(j>Y (k\q-q(s)\) u{s) ds , for q E £l\d£l . (81) 
an 

If one uses a desymmetrization, such as eq. (|72|), (|73|) or (|74|) , the above formula (|8l]) has to be 
modified accordingly. 
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Figure 20: Examples of normal derivative functions u n (s) and the corresponding 
eigenfunctions in the stadium billiard (odd-odd symmetry, a = 1.8). Here black 
corresponds to high intensity of \i> n {q)\ 2 - 
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In fig. |20]we show some examples of normal derivatives u n (s) and the corresponding eigenfunc- 
tions of the billiard, computed via eq. (|8l). The imaginary part of u n (s) is typically 5 or more 
orders of magnitude smaller than the real part. It is interesting to see that part of the structure 
of the eigenfunctions is also reflected in u n (s). For example for eigenstates with small probability 
in the region of the quarter circle also the normal derivative is small for s < ir/2. 



3.3.5 Spurious solutions I: Real Green function approach 



In certain situations and for some numerical methods it may happen that one obtains in addition 
to the true solutions of the Helmholtz equation ( |54] ) further so-called spurious solutions. This 
question is discussed in some of the papers on the boundary integral method, in particular see 
^U].|5^,|5B| and [5S,5S|. There are essentially two different situations in which they are encountered. 



The first is that one uses for the Green function instead of the Hankel function, see eq. (|6lD , just 
the real part, i.e. 



G k (q,q') = -Y (k\q-q'\) 



4 



12) 



This seems reasonable as according to (^) the Jo-Bessel function does not contribute to the 
eigenfunction. Moreover then one can work with an entirely real matrix for which the singular 
value decomposition can be computed much faster. However, when using this approach, there 
appear additional zeros (for each correct one there is one additional one) and the singular values 




k 2 100 



Figure 21: Using the real Green function (^) leads to spurious solutions (see the in- 
set) in addition to the correct eigenvalues marked by squares (compare with fig. |i7|) . 
For each true solution there is an additional spurious one (hardly visible at k 2 w 91 
and k 2 w 96). 
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Figure 22: Plot of the minima of the singular values around the eigenvalue 
k 2 = 81.93... with varying (3 using the parametrized Green function (|S3|). The 
insets show the corresponding structure of the first singular value with a logarith- 
mic vertical scale (matrix size for this computation: N = 200). 



loose their nice linear structure, see fig. [21 
parametrized Green function 



To overcome the problem of these additional zeros a 



1 



[f3J (k\q-q'\) + Y (k\q-q' 



(83) 



is used in fl58|. Thus for (3 = we obtain eq. (j82|) and for (3 = — i we get eq. (|6l|) . So using a purely 
real Green function means to vary (3 G R which changes the location of the spurious solutions but 
not those of the true ones. This is illustrated in fig. ^ around the eigenvalue k 2 = 81.93 . . . with 
(3 G [0,0.1]. Clearly on this scale the true solution does not change under variation of (3 (apart 
from the region of the avoided crossing which is due to the finite matrix size and gets smaller for 
larger N) whereas the spurious solution strongly varies with (3. For (3 = — ji with increasing real 
7 the additional zeros move away from the real axis and it seems that for (3 = — i they do not have 
any significant influence on the real axis. Still there could be cases where also for (3 = —i such a 
solution becomes relevant, but for convex geometries we have not encountered this situation. For 
an example of a non-convex geometry see section [3.3.6| . 

As an explicit example for the influence of parameterized Green function (|83|) let us consider 
the circular billiard with radius 1, where the Fredholm determinant reads (see e.g. |58|, 74]) 



D{k)= J] \-^kH[ iy {k)Ji{k) 



(84) 
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As this product converges absolutely in the whole complex £;-plane (apart from a cut along the 



negative real axis) zeros of D(k) occur when one of the factors in the product vanishes [74 



Clearly, the real zeros of D(k) correspond to the eigenvalues j m i of the circular billiard with radius 

flY 

1 and Dirichlet boundary conditions. The further zeros stem from the functions H\ (z) which 



have only zeros with Imz < f82|, and do not correspond to physical solutions of the interior 
problem. However, they can be related to resonances of the exterior scattering problem, but with 
Neumann boundary conditions |76|,[77[]. Because of the radial symmetry the S"-matrix is diagonal 
in angular momentum space 

S in = -Ell^g^ (85) 
Hl lY (k) 

and therefore the resonances are at those complex k for which 

^ {1) » = 0, (86) 



i.e. the same condition as implied by (|84|). 

If one uses the parametrized Green function ( |53"[ ) one can show (analogous to the derivation of 
eq. fl34|)) that for the circular billiard 

oo 

£><*>(*)= II [KkipJlW+Y/WMk)] . (87) 

l=— oo 

For (3 = 0, which corresponds to the real Green function (|8"2|), we get additional zeros of D^(k) 
when Y/(k) = 0. Varying (3 from zero to — i these real zeros turn complex. At first sight one might 
think that these are connected to the places with H\ (k) = 0, however numerical computations 
show that (for all studied cases) these move away from the real axis with a positive imaginary part 
and for /3 = — i one has H\ (k) = only for ImA; < 0. Thus the spurious solutions for the real 
Green function are not related to resonances of the scattering problem with Neumann boundary 
conditions. 

These examples suggest to use the full complex Green function (BTF) instead of the real variant 



82|) . Even though the numerical computation is more time-consuming for the complex case their 
advantages over choosing (^) are obvious as the variation of (3 is time-consuming as well (and 
non-trivial to implement in an automatic way). 

3.3.6 Spurious solutions II: Non— convex geometries 



Even when choosing the complex Green function (plf ) it is possible to encounter spurious solutions: 
For the circular the additional complex zeros of D(k) are sufficiently far away from the real axis, 
i.e. Im k -C so that they do not lead to problems with the application of the boundary integral 
method. However, when one considers different geometries the resonances of the corresponding 
scattering system could be closer to the real axis. This can be nicely studied for the annular 



sector billiard, see fig. 23, as the eigenvalues and eigenfunctions can be determined numerically 



with arbitrary accuracy. Using the ansatz |83], §25] 

ip(r, (f>) = [J v {kr) + cY v {kr)] sin(iA/>) {i 
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Figure 23: Boundary of the annular sector billiard for a = g7r and r± = 0.4 and 
r 2 = 0.6. 



with v = m-, m = 0,1,2,... and requiring ip(ri,<p) = and ip(r2,<f>) = gives the (implicit) 
eigenvalue equation 

J v {kr x )Y v {kr 2 ) - Y v {kr x )J v {kr 2 ) = . (89) 
For each m one gets a sequence of zeros k mn = ^E mn . 



Fig. |24| shows for the annular sector billiard with a = |^7r the first three singular values as a 



function of k 2 . The solutions of fl5S|) are marked by triangles. Clearly, there are additional minima, 
which can be associated with resonances of the dual scattering problem (for further details and 
examples of this association for the annular sector billiard see ||84j| ). In the limit of a — > tt these 
resonances are given by the eigenvalues of the circular billiard of radius r x with Neumann boundary 



conditions. For this billiard the ansatz ip(r, 0) = J m (kr) together with ^fe^ 
eigenvalue equation 



gives the 



mJ m {kr x ) - A;riJ m+1 (A;ri) = . (90) 



The circles shown in fig. |2| correspond to the solutions of (j90D and provide a very good description 
of the additional minima. 

Thus the question arises how to detect and distinguish these additional solutions. First, of 
course their existence and relevance strongly depends on the system one is studying. In many 
situations (for example convex geometries) there appear to be no complex solutions coming close 
enough to the real axis. Intuitively this seems reasonable as long as there are no trapped orbits 
outside of the billiard as these should give rise to resonances with small imaginary part. 

However, if such additional solutions exist they will show up in the 5 n statistics by an offset 
of +1 at each additional eigenvalue (unless one by chance misses the same number of 'correct' 
eigenvalues). If one has a system with such additional solutions one approach is to plot the 
corresponding normal derivative function u(s) and the eigenfunction. Usually they will behave 
quite differently for a correct eigenvalue and for a spurious solution. For example for the case of the 
annular sector billiard the normal derivative function for a spurious solution is discontinuous along 
the boundary and the corresponding eigenfunction also has contributions outside of the billiard, 
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see fig. |26|. Another test would be to use the normalization condition (|91]) for the normal derivative 
and compute the norm of the eigenfunction in the interior of the billiard. These two are the same 
for proper eigenfunctions whereas for spurious solutions they will disagree. Unfortunately, this is 
a highly inefficient method as the computation of the eigenfunction in Q is quite time-consuming. 
Instead of computing the normalization for the full billiard one could restrict to smaller subregions, 
e.g. for the annular sector billiard one could integrate over the region of the circle with radius r\ 
and check if it is different from zero indicating a spurious solution. For the annular sector billiard 
the additional zeros of the Fredholm determinant D(k) are complex as long as a < tt. Thus for 
N —>■ oo these minima will stay bounded away from zero in contrast to the minima corresponding 
to the eigenvalues. However, in practice it is not possible to check this as one has to make N too 
large to distinguish these from the correct eigenvalues. 

More generally, spurious solutions can be understood by a second look at the boundary integral 
equations. Namely, for the interior Dirichlet problem we have the single layer equation, eq. (|59|), 
and the double layer equation, eq. (|65|) . On the other hand, the single-layer equation for the outside 
scattering problem with Neumann boundary conditions at dfl is also given by the double layer 
equation (]65|) . (see e.g. |pl ,|50|1 ). As a consequence, scattering solutions of the outside scattering 
problem with Neumann boundary conditions at dQ may become relevant for real k. Namely, 
for resonances with small imaginary part they can lead to additional solutions for the double 
layer equation which are numerically indistinguishable from the correct solutions. However, these 
solutions do not correspond to solutions of the interior problem and they do not fulfill the single 
layer equation. So a possibility to distinguish spurious solutions for the interior Dirichlet problem 
is to check the validity of the single layer equation as well, which is only fulfilled simultaneously 
for correct solutions of the interior Dirichlet problem. 

A common approach (see e.g. and references therein) to incorporate this from the beginning 
is by combining the single layer and double layer equation using a linear superposition. By 
this the solutions of the outside problem with Neumann boundary conditions can be removed. 
Because of the singular kernel in the single layer equation special care has to be taken with the 
implementation. For the more difficult case of billiards with magnetic field see . 
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Figure 24: First three singular values as a function of E = k 2 of the annular sector billiard for 
a = !jj7r an d T\ = 0.4 and T2 = 0.6. The triangles correspond to the exact eigenvalues for the 
annular sector billiard, computed from eq. ( |8"5| ) and the circles correspond to the eigenvalues of 
the circular billiard with radius r% and Neumann boundary condition, determined via eq. fl90P. 
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Figure 25: First three singular values as a function of E = k 2 of the annular sector billiard for 
a = |7r and r\ = 0.4 and r2 = 0.6. The triangles correspond to the exact eigenvalues for the 
annular sector billiard, computed from eq. (|89| ) and the circles correspond to the eigenvalues of 
the circular billiard with radius r 1 and Neumann boundary condition, determined via eq. (POD 
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Figure 26: Normal derivative functions u n (s) corresponding to the correct eigenvalue 
with E = 663.88 . . . (left) and the spurious one with E = 691.77. . . (right). Here 
l\ = r\a, I2 = r%a + r 2 — 7*1 and Z3 = r\a + r 2 — r% + ar 2 . One clearly sees the 
discontinuity in u n (s) for the spurious solution. This is also reflected in the structure 
of the eigenfunction which for the spurious solution has its main contribution outside 
of the billiard. Notice that in both cases the eigenfunction has been computed 
according to (|80D inside and outside of Q. The fact that for the correct eigenfunction 
ipn{q) — (within the numerical accuracy) for q e M 2 \f2 is another test of the 
accuracy of the eigenvalue computations and eigenfunctions. 



3.3.7 Derived quantities in terms of the normal derivative function 

As the normal derivative function contains all information to determine the eigenfunction, it is 
interesting to see if this approach can be used to compute other quantities of interest. For example, 
if one wants to calculate expectation values (ip\A\ip) of some operator A in the state ip, one has 
to ensure that the eigenfunction if} is normalized, (if)\if>) = j n \ip(q)\ 2 d 2 q = 1. In principle this 

could be done by considering ((ip\if>)\ if>(q) of an unnormalized eigenfunction ip. However, an 

accurate computation of (ip\if>) using (|8~ID is quite time consuming. Fortunately, there is a simpler 
way to achieve a normalized if): If if) is a normalized eigenfunction with eigenvalue E = k 2 and 
u(s) is the corresponding normal derivative then we have the following normalization condition 
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for u(s) [p5| , [59| 



n(s)q(s) \u(s)\ 2 ds = k 2 



(91) 



If u(s) is an unnormalized normal derivative, then one obtains by 

\plk 



u(s 



§ n(s)q(s) \u(s)\ 2 ds 
an 



--u(s) 



(92) 



a normalized one. Starting with a normal derivative normalized in this way, any other quantities 
(e.g. expectation values) determined in terms of u(s) have the correct normalization. 

This is just the first example out of many highlighting the importance of the normal derivative 
for numerical computations of quantities related to eigenfunctions. For example, there are explicit 
expressions in terms of u n (s) to compute the 



normalization of if), eq. (0), |55|,|5^ 
eigenfunction ip, eq. fl8T|) 
momentum distribution 



4vrp2 



-ipq(s) 



pq(s)u n (s) ds 



(93) 



an 



and radially integrated momentum distribution |p5|.p6 

oo 

I(<p) := / ?„(r,¥>) 



r dr 



(94) 



sec 



36 for details. 



Husimi functions (see e.g. [57| , 
autocorrelation function of eigenstates p9|. 



In figs. |27|-p9|we show for the cardioid billiard examples of eigenfunctions in position space, the 
corresponding momentum distributions, the angular momentum distributions (for further details 
and examples see |p6| ) and the corresponding Husimi functions H n (s,p). The first example in 
fig. |2^ shows an example of a scarred state, i.e. an eigenstate which shows localization round an 
unstable periodic orbit p0| . Below the three-dimensional plot of the state is the corresponding 
density plot (black corresponding to high intensity) in which the localization is clearly visible. 
Also the corresponding three-dimensional plot of the momentum distribution ^567 (p) reveals en- 
hanced contributions in the directions if = tt/2, 3tt/2. This is also seen in the plot of I 567 (if) which 
shows that the probability to find the particle with momentum near 7r/2 is significantly enhanced 
compared to the mean of 1/(2tv). Another representation is the Husimi-Poincare representation 
H n (s,p) where s corresponds to the arclength coordinate along the billiard boundary and p cor- 
responds to the projection of the unit velocity vector after the reflection on the tangent in the 
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Figure 27: Three-dimensional plots of 1^567 (9) 1 2 ; 1^567 (p)\ 2 , their corresponding grey-scale pic- 
tures and the plot of the radially integrated momentum distribution I^j{^p). The momentum 
distribution |V'567(p)| 2 is concentrated around the energy shell, which is indicated by the inner 
circle. This state is localized along the shortest unstable periodic orbit, leading to an enhance- 
ment of 1^567 (i 3 ) | 2 near to if = tt/2, 3tt/2, also seen in the plot of hw^tp) near to the momentum 
direction (p = tt/2 (marked by a triangle). This localization is also clearly visible in the Husimi 
representation. 
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n — 116, odd symmetry 




ti/4 (p n/2 

Figure 28: Same as in the previous figure but for n = 116. In this case there is no 
prominent localization neither in position nor in momentum space. 
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a) n — 4042, odd symmetry b) n — 6000, odd symmetry 
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Figure 29: The eigenfunction in a) shows localization along the shortest unstable 
orbit which is also reflected in the momentum distributions and in the Husimi func- 
tion. The eigenfunction in b) is an example which appears to be quite delocalized 
both in position and in momentum space. The pictures look like those expected 
(according to the quantum ergodicity theorem) for a 'typical' eigenfunction. 
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point s. In this picture the localization around the unstable orbit is maybe most clearly seen; the 
places of high intensity are on the line p = (perpendicular reflection) and match perfectly with 
the position of the orbit. 



The second example shown in fig. 28 is an 'ergodic' state, i.e. a state which does not show 



any significant localization (as much as something like this is possible at low energies) neither in 
position nor in momentum space (apart from the localization on the energy shell). This is nicely 
reflected in the various representations. Two further examples are shown in fig. |2P| where a) is a 
higher lying scar and b) is another 'typical' state (in the sense of the quantum ergodicity theorem). 



4 Concluding remarks — or what's left ? 

There are many more issues related to scientific computing in quantum chaos which I did not 
mention in these notes. They for example include visualization techniques, programming of parallel 
computers (e.g. using PVM or MPI), or using vector computers etc. Also the more implementation 
specific aspects, including the choice of a programming language have not been discussed. A good 
starting point to learn about computing in quantum chaos are quantum maps as their numerics 
is much easier (one can use a black-box routine to get all eigenvalues at once) than for billiard 
systems, where more complicated methods have to be used. 

Acknowledgements 

I would like to thank the organizers of the summer school The Mathematical aspects of Quantum 
Chaos I in Bologna, Mirko Degli Esposti and Sandro Grafli for all their effort, especially Mirko for 
dealing with everything (including food and wine ;-) in the 'Italian way' ! Moreover, I am grateful 
to Grischa Haag for many discussions on quantum maps, and Ralf Aurich for discussions on the 
boundary element method. I would like to thank Professor Frank Steiner and Silke Furstberger 
for useful comments on the manuscript and Professors Uzy Smilansky and Andreas Knauf for 
useful remarks on an earlier version of this text. I would like to thank Fernando Perez for pointing 
out some speed improvements of the Python implementation. Most parts of this work have 
been done during my stay at the School of Mathematics, University of Bristol, and the Basic 
Research Institute in the Mathematical Sciences, Hewlett-Packard Laboratories Bristol, UK. In 
particular I would to thank Professors Jonathan Keating and Sir Michael Berry for all their 
support. Moreover, I am particularly indebted to Professor Steiner and the University of Ulm 
for their support throughout the last months. I acknowledge partial support by the Deutsche 
Forschungsgemeinschaft under contract No. DFG-Ba 1973/1-1. 



Appendix: Computing eigenvalues of quantum maps 

The first thing when thinking of solving a certain problem numerically is to decide on the pro- 
gramming language. There are numerous possibilities, ranging from Assembler, Fortran, Pascal, 
C, C++, Java, etc. to using packages like Octave, Matlab, Maple or Mathematica. Here I will 



use the quite recent scripting language Python Of course it is beyond the scope of this 

text to give an introduction to this language; several excellent introductions can be found on the 
Python homepage. In addition to the basic Python installation you will also need the Numeric 
package [|92|| , which is also simple to install. The following programs together with further infor- 
mation can be obtained from [[2(J. If you have been wondering about the name - yes it originates 
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from Monty Python's flying circus, and at several places the documentation refers to more or less 
famous Monty Python sketches. 

So here is pert_cat .py (the full version can be obtained via f20fl): 

# ! /usr/bin/env python 

import cmath 

from Numeric import zeros , Float , Complex 
from math import sin,pi,sqrt 
import LinearAlgebra 

def quantum_cat(N, kappa) : 

"""For a given N and kappa this functions returns the corresponding 
unitary matrix U of the quantized perturbed cat map. 

II II II 

mat=zeros((N,N) , Complex) # complex matrix with NxN elements 

1=1 j # predefine sqrt(-l) 

# now fill each matrix element 

# (note: this can be done much faster, see the on-line version) 
for k in range(0,N): 

for 1 in range (0,N): 

mat [k , 1] =cmath . exp (2 . 0*I*pi/N* (k*k-k*l+l*l) + \ 

I*kappa*N/2 . 0/pi*sin(2 . 0*pi/N*l) ) /sqrt (N) 

return (mat) 

def compute_evals_pcat (N, kappa) : 

For a given N and kappa this functions returns the eigenvalues 

and eigenphase of the unitary matrix U filled via quantum_cat (N, kappa) . 

H H H 

matU=quantum_cat(N, kappa) # fill matrix U_N 

evals=Linear Algebra. eigenvalues (matU) # determine eigenvalues of U_N 

# determine phase \in [0,2\pi] from the eigenvalues 
phases_N = arctan2 (evals . imag, evals . real) + pi 

# useful to determine level-spacing 

phases = concatenate ( [phases_N, [phases_N[0] +2. 0*pi]] ) 
return (evals .phases) 

### Main (used if pert_cat.py is called as script) 
if name == ' main ' : 

from string import atoi,atof 

import sys 

# Determine eigenvalues and eigenphases 

(evals .phases) =compute_evals_pcat (atoi (sys . argv [1] ) , atof (sys . argv [2] ) ) 

for k in range(0,N): # print eigenvalues 

print ("% e % e % e % e ") % \ 

(evals [k] . real , evals [k] . imag, phases [k] , abs (evals [k] ) ) 
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The only drawback of the above code is that the loop to fill the matrix is slower than a 
corresponding code in C or Fortran (notice that there are some very nice ways of overcoming this 
by inlining of code or on-the-fly compilation which are presently being developed for example in 
the context of SciPy [^3]). However, as di agonal ize uses the LAPACK library the most time- 
consuming part (at least for larger N) is done in an efficient way (not taking into account the 
possibility of using ATLAS [23] for further speed improvements). 

As a first test do (for N = 101 and k = 0.3) 



python pert_cat.py 101 0.3 



It will output the (complex) eigenvalues as a sequence x, y pairs. As a test, whether these all lie 
on the unit circle the third column is the absolute value of the eigenvectors. To plot the resulting 
data you may use 



python pert.cat.py 101 0.3 > pcat_101_0-3.dat 



which redirects the output of the program to the file pcat_101_0 . 3 . dat. To plot the resulting file 



use your favourite plotting program, e.g. for gnuplot j93f just do 



plot "pcat_101_0.3.dat" using 1:2 with points 

Now we would like to compute the level spacing distribution. To do this let us use an interactive 
Python session in which we do 

from Numeric import * # Numeric package 

from pert_cat import compute_evals_pcat # the above pert_cat routines 

from AnalyseData import * # histogram routine (see below) 



N=53 

kappa=0 . 3 

(evals, phases) =pert_cat . compute_evals_pcat(N, kappa) ; 

# sort and unfold phases 
s_phases=Numeric . sort (phases) *N/ (2 . 0*pi) 

# determine Level spacing 

# (by computing the difference of the shifted eigenphases) 
spacings=s_phases [1 : ] -s_phases [0 : N] 

(x_histogram,y_histogram)=histogram(spacings ,0.0,10.0,100) 
store_histogram(x_histogram,y_histogram, "histogram.dat") 



Then use your favourite plotting program to plot the level spacing distribution. For gnuplot you 
could do 



goe_approx(x)=pi/2 . 0*x*exp(-pi/4*x*x) 
gue_approx (x) =32/pi/pi*x*x*exp (-4/pi*x*x) 

plot "histogram.dat" w l,goe_approx(x) ,gue_approx(x) ,exp(-x) 



Here the routines to compute and store the histogram are in AnalyseData. py whose core reads 
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def histogram(data,min,max,nbins) : 
from Numeric import * 

# first select only those which lie in the interval [min,max] 
hdat=compress ( ( (data<max) * (data>min) ) , data) 

bin_width= (max-min) / nbins 

# define the bins 
bins=min+bin_width*arange (nbins) 

# determine indices 
inds=searchsorted(sort (hdat) ,bins) 
inds=concatenate( [inds , [len(hdat)] ] ) 

# return bins and normalized histogram 

return(bins , (inds [1 : ] -inds [ : -1] )/ (bin_width*len(hdat) ) ) 

def store_histogram(x_distrib,y_distrib,outdat) : 
bin_width=x_distrib [1] -x_distrib [0] 

f =open(outdat , "w") # open file for writing 

for k in range (0, len(x_distrib) ) : 

f. write (*7. e % e \n" 7, (x_distrib [k] ,y_distrib [k] ) ) 
f.writeC 10 /. e % e \n" % (x_distrib[k]+bin_width,y_distrib[k] )) 

f . closeQ 



Again, for further details and full routines see p0|. 
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